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Abstract 

We use Monte Carlo simulations to measure the spin-spin correla- 
tion function in the disordered phase of two-dimensional (;-state Potts 
models with q = 10,15, and 20 at the first-order transition point (3f 
To extract the correlation length from the exponential decay of the 
correlation function over several decades with the desired accuracy we 
make extensively use of cluster-update techniques and improved esti- 
mators. Our results for are compatible with an analytic formula. 
As a byproduct we also measure the energy moments in the disordered 
phase and find very good agreement with a recent large q expansion 
at fit- 



1 Introduction 



An important quantity to characterize the properties of a statistical system is 
the correlation length ^ which can be extracted from the exponential decay 
of a correlation function G{x) in the limit of large distances x. Usually 
various definitions of G{x) are possible and it is a priori unclear which one is 
best suited in numerical Monte Carlo simulations. There are only very few 
models for which the correlation length is exactly known and can thus serve 
as a testing ground for the employed numerical techniques. The best known 
example is the two-dimensional Ising model where the correlation length is 
exactly known at all temperatures in both the high- and low-temperature 
phase im]. Only quite recently also for two-dimensional g-state Potts models 
on simple square lattices an analytic formula for the correlation length could 
be derived p|, 0]. Here, however, the correlation length is only known at 
one special temperature, namely at the first-order transition point (3t of this 
model for q > 5. More precisely, by comparing with a large q expansion, it 
could be argued that the analytic result in Ref.||^ refers to the correlation 
length ^diPt) in the disordered phase. 

Using exact duality arguments and the (weak) assumption of complete 
wetting (which can only be proven in the limit of large q) this result was then 
converted into an explicit expression for the order- disorder interface tension, 
o'od = l/2^d 0. This formula turned out to be in good agreement with 
previous (and thus completely unbiased) numerical interface-tension data for 
g = 7 p and g = 10 [H, and also subsequent high-precision studies obtained 
compatible values |0 . The purpose of this note is to present direct numerical 
tests of the formula for ^d{Pt)- 



2 The model and observables 

The two-dimensional g-state Potts model is defined by the partition function 
i 

Z = Y.e-^^; E = -Y.5.^.^; s. = l,...,g, (1) 

{si} (ij) 

where i = {ix,iy) denote the lattice sites of a square lattice of size V = 
Lx X Ly, {ij) are nearest-neighbor pairs and 5s^s is the Kronecker delta 
symbol. In the infinite volume limit this model exhibits on simple square 
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lattices for g < 4 (g > 5) a 2nd (1st) order phase transition at /9( = ln(l + y/g). 
At pt also the internal energy densities Co and of the ordered and disordered 
phase are known exactly while for the corresponding specific heats only the 
difference Ac = — Cq could be derived analytically. 

In the disordered phase the spin-spin correlation function can be defined 

as 

Git,j)^{ds,s,-^). (2) 

For numerical purposes it is more convenient to consider the ky — projection 
of G, 

= j), (3) 

y iy,jy 

which should be free of power-like prefactors in the large- distance behaviour. 
For periodic boundary conditions translational invariance implies that g de- 
pends only on \ix — jx\, and for convenience we shall sometimes simply write 
g{x). A useful test of the consistency of our data is provided by the magnetic 
susceptibility 

x-y^^{[UQS..,i-m, (4) 



which can be computed from the area under the correlation function, 

X-y^^,M-^,M. (5) 



3 The simulation 

In our Monte Carlo study we investigated the correlation function in the 
disordered phase at Pt for Q — 10, 15 and 20 on lattices of size V — L x L 
and V = 2L X L with L = 150,60 and 40 14^^). To take advantage 
of translational invariance we used periodic boundary conditions but chose 
the lattice sizes large enough to suppress tunneling events. This guaranteed 
that, starting from a completely random configuration, the system remained 
a sufficiently long time in the disordered phase to perform statistically mean- 
ingful measurements. Since in this situation it is not obvious which update 
algorithm performs best we first performed for the L x L lattices a quite 
elaborate efficiency study of the most popular update algorithms, the local 
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Metropolis and heat-bath algorithms and the non-local Wolff single-cluster 



[^] and Swendsen-Wang multiple-cluster |jTO| algorithms. By measuring the 
integrated autocorrelation times T[^t,e of the energy it became immediately 
clear that the Metropolis algorithm is not a good candidate; see Table 1. 
Also the multiple-cluster algorithm seems to be inferior in this application. 
The other two algorithms, on the other hand, exhibit comparable rint,e; in 
particular for large q where the average cluster size is small. Taking into 
account the details of our implementation, this makes it difficult to decide 
between the two alternatives on the basis of rint.e alone. Being mainly in- 
terested in the long-distance behaviour of correlation functions we therefore 
also looked at the integrated autocorrelation times Tint,g{a;) of these quantities. 
Our results for q = 20 are shown in Fig. 1. As expected we find that for the 
local heat-bath algorithm the autocorrelations grow with distance, while for 
the non-local single-cluster algorithm they decrease. On the basis of these 
tests we finally decided to use the single-cluster algorithm for all production 
runs. It should be mentioned that for any algorithm we used the multiple- 
cluster decomposition of a given spin configuration for measurements using 
the improved estimator 

G(z,j) = ^(e(^,j)), (6) 

q 

where ^(^, j) = 1, if z and j belong to the same cluster, and G = other- 
wise. By performing the summations in eq.(0) one easily derives an improved 
estimator for g{ix,0). 

In the production runs we updated the spins after many single-cluster 
iterations with one multiple-cluster step to facilitate the most efficient use of 
the "improved estimator" (^). In units of rint,e the run time on the L x L 
{2L X L) lattices was about 35 000 (60 000) for q = 10, 116 000 (230 000) 
for q = 15, and 72 000 (35 000) for q = 20. All error bars are estimated by 
means of the jack-knife technique [|ll[]. Finally it is worth mentioning that 



all our correlation function data are stored in such a way that they can be 
reweighted to nearby temperatures in both directions; in this way we have 
also computed extrapolations of the correlation length into the metastable 
disordered region O . 



4 Results 
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4.1 Energy moments 



To convince ourselves that the system was always in the disordered phase, 
we monitored the time series of the energy measurements and computed the 
first three moments of the energy distribution, = {E)/V, q = f3fiJ,^^^ = 
(3f{{E^) - {Ef )/V, and = {{E - {E)f)/V. While can be compared 
with exact results, q and can be related by duality to the corresponding 
moments in the ordered phase. 



Cd 
(3) 



Co + /5t^(ed-eo)/v^, 

-/if + 2(1 - + 3(e, - e„)/g + Gco/A^y^, 



(7) 

(8) 



which have recently been estimated by means of Pade extrapolations of large 
q series expansions [1^ . Our Monte Carlo estimates for the Lx L and 2Lx L 
lattices can be found in Table 2, together with the Pade extrapolations as 
given in the reanalysis of Ref.|]14[ (using series expansions extended by one 
term), which are practically indistinguishable from our own Pade analysis. 
A comparison of the two sets of numbers shows excellent agreement between 
the two methods, even for the third moment and small q. 

Estimates of q from the finite-size scaling behaviour of the specific-heat 
maxima gave a consistent value of 6.0(2) [|15[ for q = 20, but much too 
small values for g = 10 |jl6|, [l^, while recent estimates from very long high- 



temperature series expansions [T^ are too large by a factor of about 2. Only 
finite-size scaling at the transition point jSt seems to give sensible results, at 
least for g = 10 [|17 . 



4.2 Susceptibility 

As a further test of the consistency of our data we compared the magnetic 
susceptibility computed according to eq.(|^) with measurements using the 
improved cluster estimators 

X={\C\)sc={\C\')sw/{\C\)sw, (9) 



where {■)sc ((■)sh/) refers to the average taken from the single (multiple) 
cluster update. As is shown in Table ^, in all cases we obtained excellent 
agreement between the three estimators. 
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4.3 Correlation function 



Let us now turn to the main subject of this note, the correlation function. A 



prehminary report of a first set of L x L data was recently given in Ref . 
Our complete set of measurements now consists of data for G along the 
coordinate axes and for the projected correlation functions g{x) for q = 
10, 15, and 20 on L x L and 2L x L lattices with L = 150, 60 and 40. The 
average of the ky = and kx = projections on the L x L lattices and the 
ky = projection on the 2L x L lattices, i.e. g{x), are shown in the semi-log 
plots of Fig. 2. The quite pronounced curvature for small x indicates that 
the simplest two-parameter Ansatz for periodic boundary conditions, 

^(x)=ach(^^/^), (10) 

which takes into account only the lowest excitation (largest correlation length), 
can only be justified for very large x. We have therefore considered also the 
more general Ansatz 

9{x) = achi^y^) + 6ch(c^^^^), (11) 

with four parameters a, 6, c, and ^d- 

Since non-linear four-parameter fits are notoriously difficult to control, 
we first fixed at its theoretical value (^^ = 10.559519..., 4.180954..., and 
2.695502... for q = 10,15, and 20, respectively), and optimized only the 
remaining three parameters. The resulting fits to the L x L and 2L x L data 
are shown in Fig. 2 as dotted and solid lines, respectively. Over a wide range 
up to about X ~ (5 . . . 6)^^ the lines are excellent interpolations of the data. 
At very large distances, however, we also see a clear tendency of the fits to 
lie systematically above the data. This already indicates that unconstrained 
fits to the Ansatz (|Tip over the same x range with as a free parameter 
should somewhat underestimate the analytical value of ^d- 

In fact, this is what we observed in the unconstrained fits to both the LxL 
and 2L X L data. In order to estimate systematic errors we performed fits to 
both Ansatze using varying fit intervals. As a general tendency we noticed a 
trend to higher values for ^d when restricting the fit interval to larger x values, 
but then also the statistical errors grow rapidly. For g = 10 this is illustrated 
in Fig. 3(a), where Xmin denotes the smallest x value included in the fits. The 
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last point used was Xmax = L/2 ioi both geometries. For the four-parameter 
fits we have stopped increasing Xmin as soon as the error on the amphtude 
h became comparable with its central value. For a reasonable range of Xmin 
values satisfying this criterion our results are collected in Table 4. Here we 
also give the results of fits of giy), i.e., the fcx = projection along the short 
direction of the 2Lx L lattices. The fits of g{x) with the smallest x^i^ values 
are shown in Fig. 2 as long [L x L) and short (2L x L) dashed lines. For 
the parameter c we obtain from the unconstrained four-parameter fits the q 
independent estimates of c ~ 1.5 — 2, with a clear tendency of decreasing 
c for increasing Xmin- This observation is consistent with the constrained 
three-parameter fits with C,d held fixed at its analytical value where we find 
the quite stable estimate of c = 1.5 ± 0.1, again for all three values of q. 

Our numerical estimates for C,d underestimate the analytical values by 
about 10 — 20% for both lattice geometries. The relative deviation clearly 
increases with increasing q. For some fit ranges we have repeated the analysis 



using so-called correlated fits pO| which, in general, seemed to be a little 
more stable. We did not observe, however, any significant increase of the 
estimates for ^d- We also investigated whether the Fourier transforms of g or 
G are less susceptible to systematic corrections and thus easier to analyze. 
Unfortunately, the answer is no. In fact, the fitted values of turn out to 
be even smaller than in the corresponding real space fits (if comparable fit 
intervals are used). 

Of course, the problem is that at the distances we have studied so far 
(a^max = L/2 ^ 7C,d) even higher excitations cannot be neglected. Due to 
convexity properties it is then natural that is underestimated by using the 
truncated Ansatz ([TT|). This general trend is illustrated in another way in 
Fig. 3(b) where we plot for g = 10 an effective correlation length defined by 
the local slopes of g{x), 

a^ = l/\n[gix)/gix + l)], (12) 

as a function of the distance x for both the L x L as well as the 2L x L data. 
For large x we expect = C,d. We do oberserve a clear increase of Q^, but 
it is of course still a long way to = 10.56. In particular with the L x L 
data is is difficult to extrapolate to the correct value since the effects of the 
periodic boundary conditions set in much too early. For the 2Lx L data, on 
the other hand, the three-parameter fit (with held fixed at its theoretical 
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value) indicates how the data should behave for very large distances in the 
long direction of the lattice. At x = 100, however, g{x)/g{0) ^ 5 x 10~^ 
which is very difficult to measure accurately, even with cluster algorithms 
and improved observables. In fact, this number reflects how improbable it is 
to generate a cluster with diameter of about 100 (recall the improved estima- 
tor (^). To cope with this problem we are currently investigating a special 
type of simulation with a reweighted Hamiltonian designed to increase these 
probabilities. Using standard simulation techniques it would take an enor- 
mous amount of computing time to follow the decay of correlation functions 
over more than 5 or 6 decades with the necessary accuracy. 

As a check of our analysis we put g = 2 in our programs, and thus 
simulated the Ising model in the disordered phase at /? = 0.71 ~ 0.80/?c- 
Here the exactly known correlation length is = 2.728865 ... [||], a value 
that is comparable to ^diPt) of the g = 20 Potts model. Our data points 
for g{x) on L X L and 2L x L lattices with L = 40 shown in Fig. 4 look 
perfectly straight in a semi-log plot. Consequently, the much simpler fit of 
the form (p!0| ) was sufficient. As a result we obtained from the fits (with 
a^min = 1) shown in Fig. 4 the estimates of = 2.7232(35) (L x L) and 
= 2.7275(24) (2L x L), and for a fit in the short direction of the 2L x L 
lattice C,d = 2.7283(20). All these estimates are in very good agreement with 
the theoretical value, showing that the employed techniques work at least in 
principle. 



5 Discussion 



Previous numerical estimates of the correlation length at (3t for g = 10 
pl| , p2| , |2^ resulted in values of about ^ ~ 6 which are much smaller than 
the theoretical prediction of ^d = 10.56. The Fernandez et al. |2^ value 
of 6.1(5) is obtained from an extrapolation of simulations aX (3 < (3t and 
thus definitely refers to the correlation length in the disordered phase. From 
our experience with correlation function fits and direct tests we believe that 
their values of ^d for (3 < [3t are already underestimated. Since the simula- 
tion points are relatively far away from /3t, the systematic errors are further 
enhanced by the extrapolation procedure used in Ref. |22|. The interpreta- 



tion of the data of Peczak and Landau I^Tl and Gupta and Irback |2^ is less 
clear to us. By repeating the simulations of Ref. we are quite convinced 
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that their technique yields a weighted average of the ordered and disordered 
correlation function which is then analyzed to obtain ^. By using a pro- 
jection to a momentum ki = 271 /L, the ordered phase is treated properly, 
but the weighted average makes the final interpretation somewhat unclear. 
Similarly, since the simulations in Ref . are performed at the specific-heat 
maximum whose finite-size scaling behaviour is governed by the transitions 
between the ordered and disordered phase, it is very unlikely that their ^ 
refers to a pure phase correlation length. In view of these problems it is 
astonishing that all three methods yield about the same value for ^. In order 
to understand this puzzling coincidence we are currently investigating also 
the correlation length in the ordered phase and first results will be available 
soon in a separate publication . 

Constrained fits with ^dif^t) held fixed at its theoretical value clearly indi- 
cate that our data for the projected correlation function (7 (x) in the disordered 
phase at Pt are compatible with the analytical prediction of Refs. 0, |] • By 
performing unconstrained fits, however, we cannot really confirm the theo- 
retical values. Rather we systematically underestimate by about 10 — 20% 
in simulations of L x L as well as 2L x L lattices. We attribute this to higher 
mass excitations which cannot be neglected at the distances investigated so 
far. To include these corrections in the fits, however, would require much 
more precise data. In a comparative study of 2D Ising correlation functions 
no such problems were encountered and the exact value of the correlation 
length could be reproduced with high precision. 
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Tables 



Table 1: Integrated autocorrelation times of the energy on L x L lattices at 
f3t in the disordered phase for different update algorithms in units of sweeps. 
The results for Wolff's single-cluster update are rescaled to these units. 



Algorithm 


q = lQ 


q = 15 


q = 20 


Metropohs 


2000 


400(100) 




Hcatbath 


125(25) 


19(5) 


11(1) 


Swendsen-Wang 


175(20) 




67(9) 


Wolff 


52(7) 


23(2) 


17(2) 



Table 2: Comparison of numerical and analytical results for energy moments 
at /3t in the disordered phase. 



Observable 


g = 10 


g = 15 


g = 20 


Crf (MC, LxL) 
Cd (MC, 2L X L) 
Cd (exact) 


-0.96812(15) 

-0.968190(81) 

-0.968203... 


-0.75053(13) 

-0.750510(65) 

-0.750492... 


-0.62648(20) 

-0.626555(97) 

-0.626529... 


Cd (MC, LxL) 
Cd (MC, 2L X L) 
Cd (large q) 


18.33(17) 
18.34(12) 
18.5(1) 


8.695(47) 
8.665(29) 
8.66(3) 


6.144(43) 
6.140(27) 
6.133(5) 


fif {MC, LxL) 
/if (MC, 2L X L) 
(large q) 


-2010(100) 

-2031(73) 

-1833(40) 


-171.0(5.1) 
-176.1(3.8) 
-174(4) 


-54.7(1.9) 
-53.9(1.5) 
-54.6(4) 
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Table 3: The magnetic susceptibility at jSt in the disordered phase from 
three different estimators. 



Lattice 


Observable 


q = lO 


g = 15 


g = 20 


Lx L 


{\C\)sc 

{msw/{\C\)sw 


38.02(14) 
38.02(14) 
38.02(14) 


10.228(19) 
10.234(19) 
10.228(19) 


5.874(11) 
5.872(11) 
5.874(11) 


2Lx L 


{\C\)sc 

{msw/{\C\)sw 


38.075(80) 
38.094(80) 
38.075(80) 


10.2330(94) 
10.2331(91) 
10.2330(94) 


5.8813(60) 
5.8808(59) 
5.8813(60) 



Table 4: Numerical estimates of the correlation length Cd{Pt) from four- 
parameter fits to the Ansatz (|TT]) in the range Xmin • • -a^max = L/2. For the 
2L X L lattices the fits along the x and y direction are distinguished by the 
index. 





g = 10 
L = 150 


g = 15 
L = 60 


g = 20 
L = 40 


Lattice 




■^min 


■^min ^d 


Lx L 


11 8.8(3) 
16 8.9(4) 
20 9.0(5) 


5 3.60(10) 
7 3.67(13) 
9 3.70(16) 


3 2.21(6) 

4 2.21(7) 

5 2.24(6) 


{2L X L), 


11 9.0(4) 
16 9.5(6) 
20 10.2(9) 


5 3.52(5) 
7 3.54(7) 
9 3.59(10) 


3 2.21(3) 

4 2.23(4) 

5 2.23(5) 


(2L X L)y 


11 8.9(4) 
16 9.1(5) 
20 9.3(7) 


5 3.52(8) 
7 3.58(11) 
9 3.62(16) 


3 2.26(4) 

4 2.30(6) 

5 2.33(7) 


exact 


10.559519... 


4.180954... 


2.695502... 
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Figure Headings 



Fig. 1: Integrated autocorrelation times of g{x) for q = 20. 

Fig. 2: Semi-log plots of the correlation functions g{x) vs distance x on 
L X L and 2L x L lattices for (a) q = 10, (b) q = 15, and (c) q = 20. 
The solid and dotted lines are three-parameter fits to the Ansatz (|l^) 
with C,d held fixed at its theoretical value. The short and long dashed 
lines show unconstrained four-parameter fits over the same x range. 
For clarity some data points are discarded. 

Fig. 3: (a) Results for of the various fits for g = 10 using all data points 
with Xmin < a; < Xmax = i^/2 as a function of Xmin- 



(b) The effective correlation length (|T2D vs distance for q = 10. The 
dashed lines are constrained three-parameter fits to the data and the 
horizontal line shows the theoretically expected result for 

Fig. 4 Semi-log plot of the correlation function g{x) of the 2D Ising model 



at (3 = 0.71 ~ 0.80/3c- The two curves are fits to the Ansatz (10) with 
= 2.7232(35) (L x L) and = 2.7275(24) (2L x L), in excellent 
agreement with the exact result Cd^"^"^^ = 2.728865 .... 
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Figure 1: Integrated autocorrelation times of g{x) for q — 20. 
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Figure 2: Semi-log plots of the correlation functions g{x) vs distance x on 
L X L and 2L x L lattices for (a) q = 10, (b) q = 15, and (c) q = 20. The 
solid and dotted lines are three-parameter fits to the Ansatz (^TJ) with 
held fixed at its theoretical value. The short and long dashed lines show 
unconstrained four-parameter fits over the same x range. For clarity some 
data points are discarded. 
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Figure 3: (a) Results for of the various fits for g = 10 using all data 
points with x-aun < a; < Xmax = as a function of Xmin- 
(b) The effective correlation length (|1^) vs distance for g = 10. The dashed 
lines are constrained three-parameter fits to the data and the horizontal line 
shows the theoretically expected result for ^d- 
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Figure 4: Semi-log plot of the correlation function g{x) of the 2D Ising 
model at /9 = 0.71 ~ O.SOPc- The two curves are fits to the Ansatz ([T0|) with 
= 2.7232(35) (L x L) and = 2.7275(24) (2L x L), in excellent agreement 
with the exact result Cd^^""^^ = 2.728865 .... 
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